ecospat包的解读

## 目前ecospat实际上提供了丰富的接口；自己实际能用到功能还是很少！

## ecospat包的开发者提供了两种脚本来模型异域分布和同域分布；、
# 参见写的那本书；

## 别的论文研究者自己开发了一种集成的思路来获取建模结果，也挺好用的；
# 如下：

########ecospat########
## 别的论文开发者提供的R代码！
library(ecospat)
library(raster)
library(dismo)
data <- read.csv('E:/Rworkspace/species_distribution/Stephanomeria.csv')                #import parental species and hybrid species distribution data (latitude, longitude and environment variables)
clim <- read.csv('E:/Rworkspace/random_distribution/NA.csv')                                 #import 10,000 random background points (latitude, longitude and environment variables)
sp.list <- levels(data[,1])
sp.nbocc <- c()
for (i in 1:length(sp.list)){sp.nbocc<-c(sp.nbocc,length(which(data[,1] == sp.list[i])))}
#Count the number of distributions for each species
sp.list <- sp.list[sp.nbocc>4]    #remove species with less than 5 distributions
nb.sp <- length(sp.list)    #the number of species
ls()
# variables selection
# try with all and then try only worldclim Variables
Xvar <- c(3:18)
nvar <- length(Xvar)
#the number of interactions between equivalence and similarity test
iterations <- 100
#the resolution of the climate space grid
R <- 100
#################################### PCA-ENVIRONMENT ##################################
all<-rbind(data[,Xvar+1],clim[,Xvar])
w <- c(rep(0,nrow(data)),rep(1,nrow(clim)))
pca.cal <- dudi.pca(all, row.w = w, center = TRUE, scale = TRUE, scannf = FALSE, nf = 2)
####### selection of species ######
sp.list
sp.combn <- combn(sp.list,2)
nsp<-ncol(sp.combn)

# storage matrices
overlap<-matrix(nrow=nsp,ncol=nsp,dimnames = list(sp.list,sp.list))          #store overlap values
equivalency<-matrix(nrow=nsp,ncol=nsp,dimnames = list(sp.list,sp.list))    #store p-values for equivalency tests
similarity<-matrix(nrow=nsp,ncol=nsp,dimnames = list(sp.list,sp.list))    #store p-values for similarity tests

# loop of niche quantification for each combination of species

for(i in 1:ncol(sp.combn)) {  
  spa<-sp.combn[1,i] #name of species a
  spb<-sp.combn[2,i] #name of species b
  clim.spa<-data[data$species==spa,Xvar+1] #env data for species a
  clim.spb<-data[data$species==spb,Xvar+1] #env data for species b

  # PCA scores
  scores.bkg<- pca.cal$li    #scores for global climate
  scores.spa<- suprow(pca.cal,clim.spa)$lisup                #scores for spa
  scores.spb<- suprow(pca.cal,clim.spb)$lisup                #scores for spb

  # calculation of occurence density
  za<- ecospat.grid.clim.dyn(scores.bkg,scores.bkg,scores.spa,100)
  zb<- ecospat.grid.clim.dyn(scores.bkg,scores.bkg,scores.spb,100)

  # overlap corrected by availabilty of background conditions
  ecospat.niche.overlap(za,zb,cor=F) 

  # test of niche equivalency and similarity
  equ<-ecospat.niche.equivalency.test(za,zb,rep=100,alternative="lower") #because it's slow, put at least 100 for final analyses
  sim<-ecospat.niche.similarity.test(za,zb,rep=100,alternative="greater",rand.type = 1) 
  # both za and zb are randomly shifted in the background (previous versions of ecospat implemented rand.type =2)

  #storage of values
  overlap[sp.combn[1,i],sp.combn[2,i]]<-ecospat.niche.overlap(za,zb,cor=T)[[1]]     #store overlap value
  equivalency[sp.combn[1,i],sp.combn[2,i]]<-equ$p.D                            #store equivalency value
  similarity[sp.combn[1,i],sp.combn[2,i]]<-sim$p.D                #store similarity 21 value

  #plot            
  #ecospat.plot.niche(za,title=spa,name.axis1="PC1",name.axis2="PC2") 
  #ecospat.plot.niche(zb,title=spb,name.axis1="PC1",name.axis2="PC2")
  ecospat.plot.niche.dyn (za, zb, quant=0.05,name.axis1="PC1",name.axis2="PC2", colz1="#CD5555", colz2="#79CDCD")
  ecospat.plot.overlap.test(equ,"D","Equivalency")
  ecospat.plot.overlap.test(sim,"D","Similarity")

  #counter
  cat(paste(i))
}

overlap             #output overlap
equivalency      #output the result of niche equivalence test
similarity          #output the result of niche similarity test